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Abstract 

We present a new technique aimed at preventing plane-wave based total 

energy and stress calculations from the effect of abrupt changes in basis set 

size. This scheme relies on the interpolation of energy as a function of the 

number of plane waves, and on a scaling hypothesis that allows to perform the 

interpolation for a unique reference volume. From a theoretical point of view, 

the new method is compared to those already proposed in the literature, and 

its more rigorous derivation is emphasized. From a practical point of view, we 

illustrate the importance of the correction on different materials (Si, BaTi03, 

and He) corresponding to different types of bonding, and to different /c-point 

samplings and cut-off energies. Then, we compare the different approaches 

for the calculation of Oq, Bq, and B' in bulk silicon. 
PACS numbers: 71.10,+x ; 71.20.Ad 



1 



Typeset using REVTjtX 



INTRODUCTION 



Since the early eighties, the Local Density Appoximation (LDA) M to Density Func- 
tional Theory (DFT) |2| has proven to be a choice tool to obtain reliable total energies 
of solids. Application of the total energy method to a solid of given (or assumed) crystal 
structure usually begins by the determination of the static equilibrium properties, like the 
lattice constant a , the bulk modulus B , and possibly its pressure derivative B . This 
determination is carried out before evaluation of other properties such as phonon spectra. 
Considerable help in this task arises from the stress theorem, derived by Nielsen and Martin 
that complements a calculation of the total energy, by giving the six components of 
the macroscopic stress tensor a aj 3, at practically no cost. 

For simplicity, we first consider a cubic structure (for which there is only one lattice 
constant). The hydrostatic pressure is related to the total energy and the volume V by : 

p = p(y\ = g H + ff 22 + ^33 = dE tot 

y ' 3 dV ' y } 

In this case, the static equilibrium properties can easily be determined either from the 
curve of energy E tot versus lattice constant a (or, in an equivalent way, versus volume V) or 
from that of pressure P versus lattice constant a. These curves can be obtained by fitting a 
polynomial to the values of E to t and P calculated for several values of the lattice constant. 

We now describe the problems that arise when a plane-wave based method is used for 
total energy and pressure calculation. 

Considering a periodic cell, the electronic wave functions should be expanded in terms 
of an infinite set of plane waves (Fourier series) at each of an infinite set of fc-points in the 
Brillouin zone. The following two approximations are introduced. Firstly, a small number of 
carefully chosen Appoints (special points) are used to sample the Brillouin zone |J . Secondly, 
the wavefunctions at each fc-point are expanded in terms of a finite basis set of plane waves 
e i(k+G)r suc ] 1 that their kinetic energy is smaller than a fixed cut-off energy E cut (Hartree 
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atomic units are used throughout the paper) : 

^\k + G\ 2 <E cut . (2) 

In principle, it is possible to approach energy convergence by augmenting the number 
of special fc-points and the cut-off energy (which corresponds to increasing the size of the 
basis set). However, working with large number of special points and plane waves requires 
a huge amount CPU time. In order to perform calculations on larger and more complex 
systems, one aims at using smaller plane-wave basis sets at each fc-point without reducing 
the accuracy of the calculation. 

Differences in the total energy of systems with the same unit cell are known to be 
accurately calculated for a number of plane waves and special fc-points smaller than those 
required to ensure convergence of the energy, provided that identical plane-wave basis sets 
are used for each calculation 0. The error due to the truncation of the Fourier series is 
systematic and cancels out. One could expect ao, Bq, and B to have similar properties : part 
of the systematic error in E tot should cancel. However, when computing energy differences 
between systems of varying size, it is impossible to use identical plane-wave basis sets due to 
the periodic boundary conditions imposed on the edge of the cell whose dimensions change 
with the volume. Instead we must choose either to use a constant number of plane waves 
Npw HI in the basis set with E cut depending upon the volume, or a constant cut-off energy 
for determining the plane wave basis set. 

The Etot curve for constant Npw is always very smooth since the set of plane waves 
scales smoothly with the lattice constant. Unfortunately, relying on this curve leads to a 
systematic underestimation of the lattice constant, as discussed in previous papers ||. On 
the other hand, constant E cut corresponds to a constant resolution in real space, which 
gives less biased results 0; but, the curve for constant E cut is ragged. One gets a set 
of disconnected micro-curves (a micro-curve being a continuous segment of E to t or P) as 
shown in Fig. 0, whereas the experimental curves are perfectly continuous. This stair-like 
variation is due to the discontinuous increase of the number of plane waves accompanying 
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the unit cell volume changes. Every time new (k + G) vectors are added to the basis set, E tot 
discontinuously decreases due to the added variational freedom in the wave function. The 
P curve also presents discontinuities for the same reasons. Moreover, at the approximate 
minimum of the E tot curve, the pressure does not vanish. 

The purpose of this paper is to analyze how these systematic discontinuities and associ- 
ated errors in the total energy and the pressure can be dealt with. In Sec. 0, we introduce the 
definitions that will be used all along the paper, including the basic concept of interpolation 
of energy. In Sec. |T|, we present accurate, but CPU time-consuming, formulas for corrections 
to energy and pressure. In Sec. |T| and |V|, approximations that make it simpler to use 



are examined. Sec. |T|is dedicated to a new technique based on a scaling hypothesis (S.H.) 



allowing to compute the energy curve for a unique reference volume, whereas Sec. [TV] treats, 



with the same notations, the corrections proposed by Froyen and Cohen |10| to pressure 
and by Francis and Payne Jll[| to energy. We emphasize the theoretical differences between 
these techniques. In Sec. [V], we present the results obtained using the scaling hypothesis for 
silicium, barium titanate, and helium. In Sec. [VT|, we compare the different techniques by 
their effect on the calculation of a , B , and B in bulk silicon. In Sec. |V1J] , we present a 
generalization to anisotropic deformations. We introduce a correction to the stress tensor 



and apply it in the case of bulk silicon. Finally, we present our conclusions in Sec. [VT I I 



I. DEFINITIONS 

Let N PW (E cut , V) be the mean number of plane waves in the basis set (see Appendix 
[A] and [0) for a given value of the cut-off energy E cut and of the volume V. The subscript 
"d" stands for "discontinuous". As the reciprocal lattice is discrete, Npw can stay constant 
for a range of cut-off energies, or for different volumes. On the other hand, it will change 
abruptly for some values of E cut , at fixed volume, or for some values of V, at fixed E cut . So, 
N PW {E cut , V) is a stair-like function (see Fig. |T|). 

Let N PW (E cut , V) be the fictitious number of plane waves in the basis set for a given 
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value of the cut-off energy E cut and of the volume V, if it were determined by the product of 

1 /2 

the density of state in the reciprocal lattice by the volume of a sphere of radius (2E cut ) 1 : 

N PW (E cuU V) = (JL\ Ivr ((2E cut ) 1 / 2 ) 3 = ^- 2 (2E cut f 2 (3) 

where the subscript "c" stands for "continuous". 

Let E c cut (N PW , V) and E^ ut (N PW , V) be the inverse of N PW (E cut , V) and 
N PW (E cut , V), respectively The first can trivially be obtained from Eq. (H) : 



2TT \ 2 /3 



1 (to'Npw 



(4) 



V 

whereas the second is defined for certain values of Npw for a given volume V, and links 



these values to a semi-opened interval [El ut , Ef 



cut I 



Let E to t [N PW , V~\ be the total energy that is calculated when the number of plane waves 
used in an actual calculation, at volume V, is N PW . 

For a given volume, this function is defined only for N PW = N PW {E cutl V). In order to 
obtain a continuous curve, we interpolate through these points, so that E to t [Npw,V~\ can 
be obtained at any value of N PW . At the present stage of the discussion, the choice of the 
interpolation scheme is not relevant (see Sec. |V|). The introduction of such an interpolating 
energy curve was also performed in other papers treating the problem [10,11]. 

We also define the function Ef ot {E cut , V} which gives the total energy obtained when the 
basis set is determined by E cut through Eq. @. This function is connected to E tot ^N PW , V] 
by: 

E d tot {E cut , V} = E tot [N PW (E cut , V),V]. (5) 

Since N PW (E cut , V) changes by abrupt jumps, this function is not the smooth desired 
one, but a set of micro-curves which correspond to a fixed Npw (see Fig. 0). 

We finally introduce an ideal continuous function E^ ot {E cut , V} which corresponds to the 
energy of an hypothetical calculation, at volume V, with N PW (E cut} V) plane waves. It is 
connected to E tot \Npw, V\ by the following equation : 

Etot {Ecut, V} = E tot [N C PW (E cut , V),V]. (6) 
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The function Ef ot {E cut , V} is the smooth function that should be used for the determination 
of material properties, whereas the one we get from an usual calculation is Ef ot {E cut , V}. 
Thus, the correction to energy aims at determining the function E% ot {E cut , V} starting from 
E?ot{Ecut,V}. 

Using the same type of notations, we can now introduce the definitions related to pres- 
sure. Let P \Npw, V] be the pressure that is calculated by using the stress theorem f3|-[5[ 
when the number of plane waves used in the calculation, at volume V, is Npw- This function 
is connected to E tot \N pw, V] by : 

'd E tol \N PW ,V} \ 
dV j_ 

/ N PW 

where the subscript Npw to the parentheses indicates that the derivative of the total energy 
versus the volume is taken at constant number of plane waves, and the subscripts V = V\ 
and N pw = N PW to the vertical bar indicate where the previous expression is evaluated. 
Let P d {E cut ,V} be the pressure that is calculated with E cut and the volume V as inputs. 
This function is connected to Ef ot {E cut , V} by : 

'dEf ot {E cutl V} 



P 



Np W ,V! 



V=Vi 

PW= n PW 



(7) 



N PW =N 1 



dE, 



cut 



V=Vi 



(8) 



In fact, this function is a set of micro-curves, each of which corresponds to fixed N pw 
(see Fig. ||). So let P c {E cut , V} be the ideal continuous pressure curve. This function is 
connected to E^ ot {E cut , V} by : 

■dEl t {E cut ,V} 



P c { E tut,Vi} 



dE, 



cut 



V=Vi 
E cut =E cut 



(9) 



The function P c {E cut , V} is the function that should be used for the determination of 
material properties, whereas the one we get from an usual calculation is P d {E cut , V}. 
Thus, the correction to pressure aims at determining the function P c {E cut , V} starting 
from P d {E cuU V}. 



6 



II. ACCURATE CORRECTIONS TO ENERGY AND PRESSURE 



The correction to energy at E cut = E\ ut and V = V\, needed to connect the output of 
a computer run Ef ot {E^ ut , Vi} with the ideal value E^ ot {E\ ut) Vi} is obtained by combining 
Eqs. (|) and (§) : 



E tot{ E Lv V i} = E tot { E Lv V i} 



+ E tot [N PW (El ut , V,) , Vi] - E tot N PW (El ut , V 1 ) , V } 



(10) 



A schematic representation of this expression is given in Fig. |^. Indeed, Eq. ( |T0"D can be 
rewritten E to t(B) = E tot (C) + {E tot (D) — E to t(E)} following the notations of the figure. The 
idea is that going from point C to B is equivalent to going from point E to D, as suggested 
by the arrows in Fig. ||. 

The correction to pressure (the correcting term is called Pulay stress [|ll|] by analogy 
with the Pulay force |13|) at E cut = E\ ut and V = V\ is obtained by deriving Eq. fllOD with 
respect to the volume at fixed cut-off energy : 

+P [N PW (E l cutl V x ), V 1 ] - P [N PW (E l cutl V X ),V X 
Npw i E lav Vi) ( dEtat \N PW V] 



V 



V=Vi 



Vx \ dN PW 

This result was obtained by using the chain rule and the following derivatives 
'dN PW {E cuU V)\ N PW (E cuU V) 

E cu t 



(11) 



from Eq. 



dV 
and 



V 



(12) 



d N PW (E cut ,V ) 
dV 



0. 



(13) 



since N PW (E cut , V) is a stair function. 

The validity of Eqs. ( |10"D and (|TT|) rests only on the choice of an interpolation scheme for 
E tot {Npw, V] ■ At this stage, in order to find Ef ot {E cut , V} or P c {E cut , V} as function of 
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the volume using Eq. ( |I0| ) or Eq. (JTJ), one should consider a few values of V, then for each 
V, choose a few basis sets, and interpolate the energy to get E tot ^N PW , V~\ . Unfortunately, 
this procedure is time-consuming. 



III. SCALING HYPOTHESIS 



We now introduce a technique that allows to make the interpolation effort only for one 
given reference volume V$. It is different from the technique proposed by Francis and Payne 
for correcting the energy flllfl , or by Froyen and Cohen for correcting the pressure [|I(J (see 
Sec. [TV]). 

We make a realistic hypothesis : for the purpose of the calculation of the correction to 
energy and pressure, the difference between energy at V and at Vq at constant E cut does not 
depend on E cut : 

Etot {Ecut, V} « E c tot {E cut , V } + f(V - V ) VE cut . (14) 

This is the mathematical expression of the principle of cancellation of errors between similar 
geometries, already mentioned in the introduction. 

By deriving this equation with respect to E cut at constant V, we get : 



dEf ot {E cut , V} 



BE, 



cut 



V=Vi 
E cu t=E\ ut 



dE^ {Ecut, V} 



dE, 



cut 



V=Vo 
E cu t=E cut 



(15) 



which means that the successive derivatives of Et ot {E cut , V} with respect to the cut-off 
energy at constant volume do not depend on V Hi~4| . This last equation can be developped 
by means of Eq. Deriving the latter with respect to E cut at constant V and using the 
chain rule, we get : 

' dE c tot {Ecut, V}\ ( dEtot [Npw, V] 



dE, 



cut 



dN 



PW 



dN PW (E cut} V) 



v_ 

N PW =N c pw (E C ut,V) 



dE, 



cut 



where the last term can be obtained from Eq. 



dN PW (E cut ,V) 



dE, 



cut 



V 

2^ 



(2E, 



cut ) 



,1/2 



(17) 



When introducing these results in Eq. (|T5|), we get : 



'dE tot [N PW ,V] 



ON 



PW 



V=Vi 
TV 



pw=N PW 



Vo ( dE tot [N PW , V] 
Vx \ ON 



PW 



V=V 
Npw=N° 



(18) 



where N PW = N PW (E] ut , V\) and N PW = N PW (E l cut , V Q ) = ^N PW . The relation ex- 



pressed in Eq. is valid for any N PW and N PW connected by : 



pw — y^ iy pw 



(19) 



Let us now consider the correction to energy. We can rewrite Eq. (pi]) in the following 



way 



E tot { E cut, Vi} 



E tot {Ecut, Vi} 
rN PW (E cut ,V 1 ) 

+ /_ 

J N d pw {E cuU V 1 ) 



dEtot Np W ,V 



dN T 



V 



V=Vi 



dN r 



N 



PW —N PW 



(20) 



Inserting Eq. (|l"8"D in this expression, we get : 

E c tot {E cuU V\} « E? ot {E cut ,Vi} 
+ 



N PW (E cuuVl ) ^ / dEtot p? pwV ] 
N PW (E cut , Vl ) Vl V 9NPW / , 



V=V 



Nr 



dN f 



(21) 



By introducing the change of variables given by Eq. fll9|), we get : 

E c tot {E cuU Vi} » Ef ot {E cut , Vi} 



+ 



^NpwiE^t^) 



dE to t[N PW ,V\ 
dNpw 



V 



V_=V 
Npw^I 



dN r 



(22) 



Finally, we have : 



Etot {E C ut, Vi} « Ef ot {E cu t,Vi} 



+E tot T^N PW (E cut ,Vi),V -Etot T^N PW (E cut ,Vi),V . (23) 
. "i J L "i 

This expression allows to correct the value of energy obtained at V\ by using the interpolating 
curve E to t \Npw, V] calculated at V . 

Let us finally consider the correction to pressure. We can get different expressions de- 
pending on the way we work starting from Eq. dID|). We can either use the scaling approxi- 
mation and then differentiate with respect to the volume, or vice versa. In the first case, we 



start from Eq. fl23|) and differentiate it with respect to the volume at fixed cut-off energy. 
This leads to : 



V d . x . ( dE tot [N PW , V] 
-^N PW (E cut , Vl ) ^ -=- 



V=V 



(24) 



Using Eq. (|18|), this can be rewritten : 

N PW (El ut , V{) 



'dE tot [N PW ,V] 



V 



V=Vi 



(25) 



N PW =N PW (El ut y 1 ) 



dN PW 

On the other side, we can start from Eq. ([□]) and use the approximation Eq. (|I~5| ) and its 
consequences. This leads to : 



la 



p 



^N PW (E l cut ,V x ),V Q 



V=V _ 



(26) 



V - c . x , ( dEtot [Npw, V] 

dN PW 



Between these different, but equivalent, expressions of the correction to pressure Eq. ([24]) is 
the most convenient because it simply needs the interpolation of E tot as a function of Npw 
for a given reference volume Vq. 



IV. COMPARISON WITH OTHER TECHNIQUES 

While the scaling hypothesis (S.H.), presented in the preceding section, starts from the 
correction to energy and proceeds by derivation to obtain the correction to pressure, Francis 
and Payne have proposed a technique for the correction to energy that starts from the 
correction to pressure proposed by Froyen and Cohen [|TU[ and proceeds by integration. Let 
us recall their results. 
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A. Froyen-Cohen correction to pressure 



The expression proposed by Froyen and Cohen for the correction to pressure is the 
following : 



2E\ ut fdE? ot {E cuU V} 



3 V\ \ dE cut j v 
The definition of an interpolating energy curve is also central in this approach. 



V=Vi 
E cu t=E\ ut 



(27) 



Let us develop Eq. (27) to compare it with the accurate correction to pressure Eq. (11 
and our approximations Eqs. fl24|) and fl26|) . By using the chain rule in Eq. we easily 
get : 



P c FC {E\ ut ,V l ) = P^El^Vx) 

KwiEUVi) fdE tot [N PW ,V] 



V=Vi 



(28) 



Vi \ dN PW 

where one term is missing compared to Eq. (]TT|). 

Let us continue the analysis. By using the scaling hypothesis in the form of Eq. ( |T5D and 
Eq. (18) respectively in Eq. (27) and Eq. (28), we get : 



P C FC {E^Vx} = P'iE^Vx} 

2El ut fdEZ ot {E cut ,V} 



3 Vx V dE cut 

pd { El cuV V l} 



V=V 

E C ut=E cut 



(29) 



V,- c ( x J dE tot \N PW) V] 
- W N PW {E cut , Vl )\ 



v 



V=Vo 

N PW =^N u r 



(30) 



,Vi) 



where clearly one term is missing in regard to Eq. (|26|). When this equation is compared 
to Eq. (|23|), we see that N PW (El ut , V\) is changed into N PW (E\ ut , V\). This is due to an 
inaccurate definition of the Pulay stress (as shown in Appendix |B|). 

Between these different expressions of the correction to pressure Eq. ( |30| ) is the most 
convenient because it just needs an interpolation of E tot as a function of Npw for a given 
reference volume Vq. 
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Before going further, we illustrate the two proposed stress correction techniques in Fig. |j. 
We draw the output of a computer run P d {E\ ut) V\} and the curves P c {E\ ut1 V{\ obtained 
respectively by applying to it either the correction given by Eq. (|24| ) or Froyen-Cohen correc- 
tion given by Eq. fl3"0|). We also draw the curve obtained by applying the following pressure 
correction : 



V 



+1" 



Vi 



Np W {E^V^Vo 



-P 



.(31) 



which consists in adding to P d {E\ ut) ]/%} the missing term in Eq. (|30|) in regard to Eq. ( |26|) 
(which is equivalent to Eq. (0)). The graph clearly shows the importance of each term of the 
proposed corrections. The correction of Froyen and Cohen given by Eq. (^) is responsible 
of a shift of the uncorrected curve. Whereas, the correction given by Eq. (|31|) is responsible 
of the cancellation of the jumps between the micro-curves. The S.H. correction includes 
these two effects as it can be seen from Eq. (EjJ) and Fig. |j. The first is definitely the most 
important, but the second is not negligible. Moreover, Eq. ([24]) or Eq. fl3"0|) are equally easy 
to use, so that the more accurate Eq. (El) should always be preferred. 



B. Francis- Payne correction to energy 

The expression proposed by Francis and Payne for the correction to energy is the following 
(see Appendix |C]): 

E tot,FP { E luU Vi} = E tot { E lut, v i) 

-^k\n( = PW {ELt > Vl) ^ ( 9EL {EcutV} 



N PW (El ut , Vy) y V dE, 



■'cut 



V=Vi 
E C ut=E cut 

(32) 

It can be shown (see Appendix 0) that the derivation of this expression implies the use 
of the approximation given by Eq. dl5|) . Moreover, it is obtained by integrating the Pulay 
stress expression proposed by Froyen and Cohen which is inaccurate (see Appendix [BD . By 
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using the chain rule in Eq. 



we get : 



E tot,FP { E luti Vi) = E tot{ E Lt V l} 

-N PW (El ut , V,) In 



N PW (Ej t , V± 
N C PW V,) 



x 



'dE tot [N PW ,V] 



dN 



PW 



V=Vx 

Npw=Np W (ElutVi) 



By using our approximation in the form of Eq. (|T|) in Eq. (|32|) , we get : 



E tot,Fp{ E Lt' V l} = E tot{ E luV V l} 

-1 



2E cut j n / N PW ( E luV V l 



N PW (EU V x 




dE? ot {E cut , V} 



dE, 



cut 



By using our approximation in the form of Eq. ( PD in Eq. ( |55| ) , we get : 



v 



(33) 



V=V 
F C ut=Fl 



(34) 



E tot,FP { E luf V l} = E tot { E luti Vl] 



ttN pw (E cut , V x ) In == — 



'<9£ fot [n pw ,v] 



dN 



PW 



V=V 



(35) 



This last expression of the correction to energy is the most convenient because it just needs 
an interpolation of E to t as a function of N P w for a given reference volume Vq. 

The comparison with S.H. technique is not as straightforward as for correction to pres- 
sure. But it is clear that the inaccuracy in the definition of the Pulay stress, from which 
Eq. (p5|) is derived, is a source of error. 



V. APPLICATIONS 

The importance of the correction is now investigated for silicon, barium titanate, and 
helium, various materials presenting different types of bonding. 
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Our calculations, performed within the local density approximation (LDA) |T|, use a pre- 



conditioned conjugate gradient algorithm [15, To]. We use a rational polynomial parametriza- 



tion of the exchange-correlation energy functional |T7|], which is based on the Ceperley-Alder 
gas data [IJ|. The Brillouin zone is sampled with different Monkhorst-Pack || meshes of 
special fc-points. The "all-electron" potentials are replaced by ab initio, separable, norm- 
conserving pseudopotentials, as described below. 

Then we apply Eqs. ( p3|) and ( 24|) to correct the results. The interpolating scheme 



for obtaining E to t \Npw, Vo\ is the following. We first calculate the total energy for three 
different values of Npw '■ 

N PW (E cut - 3%, V ) , Np W (E cut , V Q ) , N PW (E cut + 3%, V ) . (36) 

Then we interpolate between these values by an exponential fit of the form : 

Etot [N PW , V ] = EZl + exp (a + aJj PW ) . (37) 

First, we apply this correction to bulk silicon, with covalent bonding. The corresponding 
pseudopotential was built following the scheme proposed in Ref. (TP]]. The atomic positions 
in the unit cell are completely determined by symmetry : the only free parameter is the 
lattice parameter a. We compute simultaneously the total energy E tot (a) and the pressure 
P{a) as we vary the lattice parameter a at a constant energy cut-off (E cut =Q Ha with 2 
special fc-points in the irreducible Brillouin zone). Then we apply S.H. correction. The 
effect of correction to energy, in Fig. [5] is obvious : it cancels out the jumps between micro- 
curves. The effect of correction to pressure in Fig. [| has already been commented. Firstly, 
the jumps between micro-curves are suppressed. Secondly, the whole curve is shifted towards 
larger pressure. These results will be commented with more details in Sec. 0. The lattice 
parameter value is found to be 10.23 Bohr. 

Then, we consider barium titanate (BaTiOa), for which we used extended norm- 
conserving pseudopotentials given in Ref. ||20|| . This material is sometimes classified as 
an ionic compound. As pointed out recently, it also presents a partial covalent character 
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21| , |22| . Its structure is cubic perovskite at high temperature (above about 120°C) while 



it undergoes 3 successive ferroelectric phase transitions as the temperature goes down. In 
this study, we will focus on the cubic, high symmetry phase in which the atoms are at 
symmetric positions so that the only structural degree of freedom is the lattice parameter 
a. When using the S.H. correction, a first estimation of a at 7.53 Bohr is already possible 
on the basis of only 4 points at 20 Ha cutoff with 4 points in the irreducible Brillouin zone 
(see Fig. ||). For BaTiC>3, an accurate investigation of some properties (as, for example, 
the phonons frequencies) requires nevertheless to work at a 45 hartrees energy cut-off on a 
6x6x6 mesh of fc-points. At this high cut-off, the total energy is well converged ( ®J? tot is 
very small) so that the correction on E tot becomes negligible (Fig. [j]). We predict a lattice 
parameter of 7.45 Bohr in good agreement with other LDA calculations. By contrast, even 
for this case, the error on the pressure remains critical. As illustrated on Fig. |7|, it can be 
efficiently corrected when using the S.H. technique. This correct estimation of the pressure 
reveals essential for the investigation of the energy surface of the low symmetry phases. 

Finally, we perform an energetic calculation for FCC helium. Once again, we observe the 
cancellation of the scattering in the total energy data and the shift of the pressure curve, 
showing that the correction is independent of the type of bonding, and applicable to a large 
range of materials. These results will be published in a future paper |[23||. 



VI. COMPARISON OF THE TECHNIQUES 



It has been shown in the previous section that corrections to energy and stresses can be 
rather important. Now, we analyze in more detail their importance, and we compare more 
quantitatively the two possible correction techniques in the case of silicon (Eqs. (|23|) and 
([24]) for the scaling hypothesis technique, Eqs. ([M]) and ( |35|) for Froyen-Cohen technique and 
Francis-Payne technique). We investigate different aspects of this correction, more specif- 
ically we analyze its effect on the calculation of the lattice constant a , the bulk modulus 
B , and its pressure derivative B . These properties can easily be determined either from 
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the curve of energy E to t versus lattice constant a (or, in an equivalent way, versus volume 
as V oc a 3 ) or from that of pressure P versus lattice constant a. 

In order to obtain the E to t{a) and P(a) curves, we calculate the total energy and the 
pressure for a set of lattice constants at a given cut-off energy. 

As we mentioned above, the result is a set of micro-curves : Ef ot {E cut , a} and 
P d {E cut , a}. By fitting a polynomial to these values, we get continuous curves Et A(i) 
and P(a). Depending on its degree, the polynomial fit matches more or less accurately the 
data. We measure the matching by the standard deviation x of the data (xi,yi) from the 
polynomial g(x) : 



X 



N ■ ■ x2 



s (yi-g( x i)Y 

N *-i (38) 



where N is the number of data. For example, from Fig. we have calculated the value of 
X for respectively uncorrected data (open circles) and those obtained with S.H. correction 
(solid circles), with respect to their corresponding polynomial fit (respectively, the broken 
curve and the solid one). We illustrate in Fig. ||the evolution of x as a function of the degree 
of the polynomial. On one side, it is evident that the higher the degree of the polynomial is, 
the better the matching with the data will be. On the other side, the aim of the fitting is to 



get rid off of the noise |24] due to jumps between micro-curves. The higher the degree, the 
greater the part of this noise included in the fit. So, a compromise has to be reached. We will 
consider that the separation between trustable data and the noise has been accomplished 
as soon as the standard deviation reaches a plateau (called residual noise) as a function of 
the degree of the polynomial. Going to higher-order polynomial would mean beginning to 
include the residual noise in the fit. So, the degree of the polynomial is chosen by detection 
of a plateau in the standard deviation of the data from the polynomial. In Fig. §, this 
plateau is already reached for a polynomial of the third degree. 

We apply the different techniques presented in the previous sections to correct the cal- 
culated values. By fitting a polynomial to these values, we get two different expressions 
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of Ef ot {E cut , a} and of P c {E cut , a}, corresponding to the Froyen- Cohen- Francis-Payne cor- 
rection or to the scaling hypothesis. So, finally, we have three different expressions of the 
E to t(a) and P(a) curves. 

However, if there are not enough data, statistical fluctuations occur on the static equi- 
librium properties (ao, -B , and B' ) due to the noise [23|| present when the correction is not 
used (see Figs. ||-[7| of Sec. |V| especially for BaTiOa). So, from now on, we will work with a 
large number of data points in order to suppress these fluctuations, and concentrate on the 
techniques. 



A. Scale Analysis 

In order to analyze quantitatively the reduction of noise, we consider different sets of 
lattice parameter, each set being characterized by a given scale which is defined as the 
distance between two successive points of the set. We consider six different scales : 0.01, 
0.025, 0.05, 0.1, 0.2, and 0.5 Bohr (see Table |). We generate the three different E tot (a) 
curves and of the three P(a) ones for all these sets. From these curves, we determine the 
lattice constant a , the bulk modulus B and its pressure derivative B . 

Many observations occur. The residual noise in the energy and pressure data is decreased 
by the correction. This noise reduction is more important for small scales than large ones 
(see Fig. ||). This is due to the fact that the presence of micro-curves is much more apparent 
at small than at large scales. So, for the energy curve, the results obtained at large scales 
with and without correction do not differ significantly. Whereas, for the pressure curve, the 
correction also includes a shift of the uncorrected curve : whatever the scale the correction is 
always important. So, from now on, we just discuss the results obtained with the corrected 
values of energy and pressure. Regarding the lattice constant, it can be seen (see Table |T[) 
that there is a quite good agreement between the values obtained at the different scales and 
between those obtained starting from the E to t(a) curve and those calculated from the P(a) 
one. Regarding the bulk modulus, there is still a very good agreement between the values 
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obtained at the different scales from the pressure curve but not between those calculated 
from the energy one (see Table [TTTD . Moreover, the values obtained starting from the E tot (a) 
curve and those calculated from the P{a) curve agree only at large scale. 

Regarding the derivative of the bulk modulus versus pressure, the agreement is worse 
(these data are not reproduced here). 

As a general rule for the calculation of static equilibrium properties (ao, -Bo, and B' ), we 
can state that firstly, calculations based on the pressure data are more accurate than those 
based on energy data, and secondly, the accuracy decreases with the number of derivatives 
taken from the starting curve. 

Regarding the different proposed techniques, it can be noticed that their results are 
nearly the same. This is due to the fitting operation that eliminates the residual noise. Note, 
nevertheless, that it is always lower with the S.H. technique, especially for the pressure curve 
(see Fig. |). 

B. Micro-curve Analysis 

When studying the small scales, in the previous section, micro-curves corresponding to 
a constant number of plane waves have been detected. The analysis of these micro-curves 
separately will emphasize the efficiency of the different corrections. 

For this purpose, we consider the scale 0.01 Bohr with the following sets of points for 
each micro-curve (see open circles in Fig. |5|) : from 10.00 to 10.04, from 10.05 to 10.12, from 
10.13 to 10.21, from 10.22 to 10.29, from 10.30 to 10.37, from 10.38 to 10.44, and from 10.45 
to 10.50. We generate for each of these micro-curves the three different E tot (a) curves and 
of the three P{a) ones. From these curves, we also determine the lattice constant ao, the 
bulk modulus B , and its pressure derivative B . 

Regarding the lattice constant, the agreement between the values obtained for the dif- 
ferent micro-curves (see mean value fi and standard deviation o in Table [IV]) is much better 
with S.H. technique than with Froyen-Cohen or Francis-Payne technique. Note also that 
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there is a very good consistency between the values obtained starting from the E to t{a) curve 
and those calculated from the P(a) one. 

Regarding the bulk modulus, there is still a very good accordance between the values 
obtained for the different micro-curves from the pressure curve but it is a little bit worse for 
those calculated from the energy one (see mean value \i and standard deviation a in Table |V|). 
Moreover, the values obtained starting from the E to t{a) curve and those calculated from the 
P{a) one agree very well. 

Regarding the derivative of the bulk modulus versus pressure, the agreement is worse, ex- 
cept for results obtained by calculations based on the pressure curve with the S.H. technique 
(these data are not reproduced here). 

The conclusions, regarding the calculation of static equilibrium properties (a , B , and 
-Bq), drawn in the scale analysis still apply for micro-curve analysis. 

Regarding the different techniques, the results are generally not the same. This is due 
to the fact these techniques transform the original micro-curve into different micro-curves. 
The agreement of the results between the micro-curves is a mesure of the effectiveness of the 
correction method. In this case, it clearly shows that the S.H. is the best. When looking at 
Fig. |, we see that there are still discontinuities (this should be analysed by further studies) 
in E to t{a) curve though there is a good agreement of the slopes between the micro-curves. 
This means that the micro-curves are parallel, which is confirmed by the fact that there are 
no discontinuities in the P{a) curve. 

We have also tested the validity of the approximation given by Eq. flT5l). The first test 
consists in using a reference volume different for each micro-curve. Doing so, the differ- 
ence between V and Vq is smaller and thus the influence of the approximation should be 
reduced. It appears that the results obtained (not reproduced here) are almost not affected, 
in favor of the scaling hypothesis Eq. (|T5|). The second test consists in not introducing the 
approximation. This can only be done for a reduced set of lattice constants, for example a 
micro-curve, because it necessitates to calculate E tot pVpw, V\ for each of the lattice con- 
stants. Unfortunately, this increases the residual noise so that no comparison can be made, 
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and this second test is inconclusive. 



C. Convergence Analysis 

We finally analyze the effect of the correction on the convergence of the calculated values 
of do, Bq, and B in cut-off energy (working successively at 3, 6, 10, and 15 Ha, the latter 
considered as giving completely converged values) and number of special fc-points (working 
successively with 2, 6, and 10 special points, the latter considered as giving completely 
converged values) for three different scales (0.01, 0.1, and 0.5). We see that working at 
E cut =10 Ha with 2 special fc-points can be considered as enough converged (it leads to a 
relative error of 0.1% on do, 1% on B , and B ). From Tables and |VI| we see that the 



convergence is not really improved by the correction except for the result obtained from 
the pressure curve. This can be explained by the fact that the fitting operation acts as a 
correction by eliminating the residual noise on the energy curve, so that the results with and 
without correction are nearly the same. In the limit of a very large number of data points, 
this implicit correction is sufficient. However, in practical applications, very few points are 
used, and the correction is needed to avoid the statistical fluctuations mentioned above. 
In this sense, it can be said to improve the convergence. For the pressure, the correction 
also includes a shift of the uncorrected curve, this explains why the correction is always 
important whatever the degree of convergence. 



VII. ANISOTROPIC DEFORMATIONS 

We now consider anisotropic deformations and generalize our technique for stress correc- 
tion. In this case, the total energy does not only depend on the volume of the unit cell but 
also on its shape. In order to be as general as possible, we consider that E tot depends on 
the matrix A formed by the components of unit cell vectors. So the energy correction given 
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(39) 



by Eq. fliCf ) becomes in the case of isotropic deformations : 

E tot { E lut > ^ } = ( E L > A 1 1 

+£ tot pV^ (^ ut , ^i),^ 1 ] " 3w [iVp W Vi),^ 1 

where V\ = detA 1 is the volume of the unit cell. We now make the approximation that the 
difference between energy at A and at A? at constant E cut does not depend on E cut : 

Et t {Ecut, A} « E c tot {E cuU £} + f(A - 4 ) V£ cui (40) 

which is a generalization of Eq. flU} ). This leads to the generalized forms of Eq. (|15D : 



'd Ej ot {E cut ,A} ' 
dE cut 

and of Eq. ([TSj) : 

'dE tot [Np W ,A\ 



A=A X 

E cu t=El ut 



'dE° tot {E cuh A} 



dE, 



cut 



A=A° 

E cut =E\ ut 



(41) 



dN 



PW 



Npw=N PW 



Vo ( dE tot [N PW ,A\ 



PW 



A=A° 

N PW =N PW 



(42) 



where 



_o V —i det^ ,,! 
* PW ~ V, PW ~ detA 1 PW ' 

Finally, the generalized energy correction writes : 
Etot {Ecut, A}} « Et{E cut) A 1 } 



(43) 



+E-, 



tot 



T7~N PW (Ecut, Vi),A° 
Vi 



E 



tot 



^N PW (E^VJ^A 



• (44) 



Deriving this expression with respect to the strain e Q/ g at fixed cut-off energy, we get the 
stress correction : 



) Vo N d (F 1 ,{ dE t ot[N PW} A] 

+Oaf3TT2^PW {E cut , Vl) 



1 V oN PW 

where we have used the following definitions : 

'dE? ot {E cut ,A}\ 

' Ecut 



A=A° 

Npw^Np^E^V,) 



(45) 



rip A 1 } 



l=o 
A^A 1 

Ecut=E cut 



(46) 
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and 

Ecu 

We have also used the following result : 
1 dV 



(47) 



E cut =El ut 



e=0 = $ a p • (48) 
V=Vi 



V 1 dt aP 

which is obtained by deriving the expression of the volume V of a unit cell defined by matrix 
A in terms of the strain tensor e with respect to the unit cell defined by matrix A 1 : 

V = detA_ = Vi det{5 + 1) = ^Vie ijk ei mn (5 u + e a )(S jm + e jm )(5 kn + e kn ) . (49) 

The third-rank tensor e is defined in order that Eij k is equal to +1 if {i,j,k} corresponds 
to any even permutation of {1,2,3}, to -1 if {i,j,k} corresponds to any odd permutation of 
{1,2,3}, to in every other case. 

It can be shown from Eq. ([15]) that the stress correction only concerns the diagonal part 
of the stress tensor a a p. This result had already been mentioned by Vanderbilt [^5J, but 
in a form similar to that proposed by Froyen and Cohen [pTOfl - More precisely, the stress 
tensor a a p can be decomposed into an isotropic contribution cr^ (where the h superscript 
stands for "hydrostatic") and an anisotropic one cr^g (where the d superscript stands for 
"deviatory") : 



a 



h 

a/3 



-PS Ql 3 (50) 



and 



a % = <*ap ~ ( 51 ) 

where P is the pressure defined by Eq. ([I]). From this decomposition, it appears that only 
the isotropic part a^p of the stress tensor has to be corrected. So, we come back to the case 
of isotropic deformations and the correction to pressure is given by Eq. ( p4| ) . 

We apply this correction in the case of silicon. The chosen anisotropic deformation 
corresponds to a compression (or expansion) of the unit cell along the [001] direction, with 
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concurrent expansion (or compression) of the unit cell along the [100] and [010] directions, 
such that the length of the three cubic directions are changed to a = b and c, while the 
volume V is unchanged. The energy and the pressure are functions of the ratio y = b/c. 
&n = 022 and (J33 stresses are present, while non-diagonal stresses vanish. Note that the 
atomic positions in the unit cell are no longer completely determined by symmetry. The 
results obtained are shown in Fig. [l(| The minimum of the total energy and the zero of 
stresses nearly correspond to the expected value of y = 1 (the relative error is 0.3% for the 
value obtained from the energy curve and 0.03% for the values obtained from stresses curves). 
The effect of energy correction in Fig. |H] is obvious. In the stress graph, first, there is a 
suppression of jumps between micro-curves. Though it is not visible on the graph, it can be 
detected by the residual noise reduction (it is approximatively divided by 3). Second, there 
is a shift of the uncorrected curve. It is the only effect included in Froyen-Cohen technique. 
In this particular case, it is a shift by a constant due to the fact that the deformation is at 
constant volume. Indeed, N c pw (E cut , V) is constant in Eq. ( pOf ) for Froyen-Cohen technique 
and in Eq. (f2|) for S.H. technique. 

VIII. CONCLUSIONS 

In this paper, a new means of correcting total energy and stresses calculations performed 
using a fixed cut-off energy plane wave basis set has been presented. This technique relies 
on the interpolation of the energy as a function of the number of plane wave, and a scaling 
hypothesis (S.H.) that allows to work with a unique reference volume. It has been compared 
to that of Froyen and Cohen for stress correction and that of Francis and Payne for energy 
correction, both presented using the same notations. On a theoretical point of view, the 
approximations used in the different approaches have been contrasted, showing that the 
scaling hypothesis technique is more rigorous. On a practical point of view, we have shown 
the importance of the correction by presenting its effects on different materials (Si, BaTi03, 
and He) corresponding to different types of bonding. 
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Then, we have compared the different methods by analyzing their effect on the calculation 
of do, -Bo? an d B' in bulk silicium starting from both the energy curve and the pressure one. 
As a general rule, calculations based on the pressure data are more accurate than those based 
on the energy data. Also, the accuracy reduces with the number of derivatives taken from 
the starting curve (i.e. pressure derivative of the bulk modulus is more difficult to correct 
than the bulk modulus). The results of a scale analysis shows that if we do not focus on one 
micro-curve, the benefits of the different correction techniques are rather similiar. This due 
to the final fitting operation that cancels the residual noise, although the latter is always 
obtained lower with the S.H. technique. By contrast, the results obtained when working at 
the micro-curve level show that the S.H. is the best, though there are still discontinuities in 
E tot (a) curve. 

Regarding the convergence with respect to the cut-off energy and to the number of special 
fc-points, we see that it is not really improved by the corrections except for the results 
obtained by the pressure curve. For the energy curve, in the limit of a very large number of 
data points, the fitting operation acts as a correction by eliminating the residual noise, and 
the results with and without correction are nearly the same. However, when working with 
a small number of data, the correction is needed to avoid statistical fluctuations. For the 
stress curve, the correction also includes a shift of the uncorrected curve, this explains why 
the correction is always important whatever the degree of convergence. 
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Since the number of plane waves differs at each fc-point, while the techniques exposed in 
this paper use a single number to characterize this set, we decide to work with an average 
number of plane waves. It is not clear how this average should be calculated on the different 
Appoints. We can think to an arithmetic mean : 



APPENDIX A: MEAN NUMBER OF PLANE WAVES 




(Al) 



k 



or to a geometric mean (as suggested by Francis and Payne [TT|) : 



N PW (E mt , V) = J] (N d PWtk (E cut , V)) 



fk 



(A2) 



k 

where are weights defined in order that : 




(A3) 



k 
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An interesting starting point for our reflexion would be the expression of the total energy 
as function of its value at the different Appoints. Such an expression is known for some of 
its components. For example, the kinetic energy is the arithmetic mean of its value Eki n ,k 
at the different /c-points : 

Ekin = fkEkin,k (A4) 

k 

where are the weights of the different fc-points. On the other hand, the exchange- 
correlation contribution to total energy is global; it can not be physically divided in a sum 
of contributions of each fc-point. However, we can suppose that Eq. (|A4j) can be generalized 
for total energy, so that we can write : 



Ef ot {E cuU V} = J2 fkKt,k {Ecut, V} , (A5) 

k 

By using Eq. @ to develop Eq. (|A5|), we get : 

E to t 



N PW (E cut , V),V\=Y1 fk E tot,k [N d PW ,k (Ecut, V),V}. (A6) 

k 

This expression can be considered as definition of the average number of plane waves. 

Unfortunately, it is not very useful as E tot ^ is not known as a function of Npy/ at each 
different fc-point. To say something about the definition of Np\y, we have to consider cases 
where the total energy at each different /c-point and the total energy E to t are the same 
function of N PW (it corresponds to deleting the k subscript of E tot in the right-hand side 
of Eq. (TA"6|)). This precisely the case where the bands are dispersionless, which is the case 
e.g. when studying an isolated molecule in a supercell. Let us analyze two of these simple 
cases. If all these functions were the same linear function of Npw, the definition of N py/ 
given by Eq. (|A6|) would correspond to an arithmetic mean. Whereas, if all these functions 
were the same logarithmic function of Npw, the definition of Npw given by Eq. (|A6|) would 
correspond to a geometric mean. 
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In our calculations, we decide to use the geometric mean for defining Np\y- This choice 
seems more physical than the arithmetic mean, because the latter suppose a linear depen- 
dence of E tot versus N PW . But it is clear that the problem should be investigated further, 
and that it is a source of error on the corrected values. 



APPENDIX B: FROYEN-COHEN TECHNIQUE FOR PRESSURE CORRECTION 



The starting point of this technique is the definition of the Pulay Stress. It is the quantity 
to be added to the pressure that is obtained using E cut and V for the definition of the plane- 
wave basis set, to get the pressure that effectively corresponds to these given values of E cut 
and V : 



a Pulay (El t , V) = P c {El t , V}-P d {El ut , V l ) 

= P c {El ut , V 1 }-P \n pw (E l cutl V X ),V X 



Froyen and Cohen, in their paper ||10|| , propose : 

2E] ut ( dEt ot {E cut ,V} 



O Pulay, FC [E^, Vl) 



dE c . 



V=Vi 



(Bl) 



(B2) 



3 Vi \ uni cu t / \ 
We now show that this expression corresponds to an inaccurate definition the Pulay stress. 
Using the definition of E c cut [Npw, V) given by Eq. (f|), we get that : 



d v , -- 3 (^™)'(vr 3 

' N PW 

The Pulay stress, as given by Eq.( P2] ) can thus be rewritten as : 

' dE c tot {E cut ,VY 



V 



(B3) 



VPulay,FC {E^t, V) 



dE, 



cut 



V 



X 



'dE c cut (N PW ,V) 



dV 



v=v Y 



N PW 



V=Vi 

N PW =N c PW (El ut ,V 1 ) 



(B4) 



Finally, the expression of the Pulay stress given by Eq. ( |B4[ ) can be worked out, using the 
chain rule and the definitions Eq. (|7|) and (|9|) of P \N P w, V] and P c {E cut , V}, in order to 
give : 

VPulay,FC Vi) = P° Vi} - P [Kw (#L, ^i), Vj . (B5) 
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This expression does not correspond to the accurate definition of the Pulay stress given by 
Eq. (|ET|), since N PW (E cut , V) has been replaced by N PW (E cut , V) in the last term of the 
right side of this expression. 



APPENDIX C: FRANCIS-PAYNE TECHNIQUE FOR ENERGY CORRECTION 



The technique proposed by Francis and Payne consists in integrating the Pulay stress in 
the form of Eq. (|B5|) and Eq. ( |B2"D between V 2 and V\ defined on Fig. [|, i.e. by the following 
relations : let V± be the volume at which energy has to be calculated, N PW = N PW (El ut , V\), 
N PW = N PW (El ut , Vi) and V 2 be the volume such that N PW = N PW (El ut , V 2 ). 

Let us now consider each terms of this integration. The first term of the right side of 
Eq. ( p5|) is quite easy to integrate : 



/ P c {El ut ,V'} dv> = - 

JVo JVo. V 



QE- ot {E cut y} 

dV 



y 2 v 7 E cu t 

E m \ E cuti ^2} - E c tot {E mt , Vi} 



V=V dV 
E C ut=E cut 



E; 



tot 



N 



PWi ^2 



E, 



lot 



N PW , Vx 



(CI) 



It corresponds to E tot (A) — E tat (B) on Fig. |3|. 

The integration of the right side of Eq. (|B2|) necessitates the use of the approximation 
given by Eq. (|l^) : 



f Vl ^ E cut ( 9E- ot {E cut ,V} \ 

Jv 2 * V ' V 9V 



v=v 



dV 



2El utln fV 2 \ fdE? ot {E cut ,V} 



3 V.V V dE cut 
Using the definition of V\ and V 2 with Eq. @, we get : 



V=Vi 
Ec U t=E\ ut 



(C2) 
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Npw (EU V 2 ) 
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N PW i E luti V h 




(C3) 



v n pw ( E Lt >Vi) J \ n pw ( E Lt ) Vi , 
A new approximation is needed to integrate of the second term of the right side of 

Eq. (Pop, because it is impossible to solve analytically : 

V 2 JV 2 V / N f 



v_=v _ 

N Pw=N c PW (El ut y' 



dV 



(C4) 
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Indeed, the derivative in this integral has to be taken on a curve at constant number of 
plane waves Npw = ~N°pw ( E cup V)i the latter varies while integrating (because of the 
variation of V. Graphically, it means that this derivative has to be evaluated on each of 
the curves at constant N pw situated between the points A and B on Fig. [3| 

On the other side, if we accept that the number of plane waves is kept constant to 
N PW = N PW (El ut , V 2 ) = N PW {El utl Vi) = N PW while integrating, we get quite easily : 



v 2 



P\N PW X] dv> =E tot 



N PW , V 2 



-E-, 



tot 



NpwM 



(C5) 



It corresponds to E tot (A) — E tot (C) on Fig. [| It should be noted that making this approx- 
imation leads to a result equivalent to the one that would be obtained by integrating the 
Pulay stress in the form of Eq. (|6"T| ) (instead of Eq. ( p35| ) ) and Eq. ( |B2| ) . 

Globally, the integration of the right side of Eq. (|B^ ) gives E tot (B) — E tot (C), which is 
the correction for energy : 



E 



lol 



N PW , V 1 



E 



lot 



N PW , V x 



.^k\xx 



N 



PW 



' dE c tot {E cuU V} 



N 



pw , 



BE, 



cut 



V 



V=Vi 



(C6) 



This can be rewritten as : 



E tot { E \ut > V i } = E tt { E lut » V i } 



2E Lt , ( N PW ( E luti V l) \ ( ® E tot { E cut, V} 



N PW (E l cut , V, 



BE, 



cut 



V=Vi 



(C7) 



Ecut-El ut 



which is the final result presented by Francis and Payne. 
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FIGURES 

FIG. 1. Number of plane waves N PW in a basis set defined by Eq. (|2|) in a 2D case. The upper 
graph illustrates how these k-points are enclosed in a circle of radius (2E cu t) l l 2 . The lower graph 
illustrates the stair-like evolution of N PW as a function of E cut . 

FIG. 2. Total energy (top) and pressure (bottom) of Si calculated at constant cut-off energy 
(E cu t=6 Ha with 2 special /c-points in the irreducible Brillouin zone [IBZ]). The graph highlights 
the presence of "micro-curves" (see text). 

FIG. 3. Schematic representation of energy correction. The desired curve Ef ot {E^, ut , V} is 
the one going through points A and B. The calculated curve Ef ot {E^ ut , V} is the one going 
through points A and C. The correction at volume V\ moves point C, obtained by a calcu- 
lation with N PW = N PW (El ut , V\) , to point B, obtained by an hypothetic calculation with 
N PW = N° PW (E^ ut , V\\ . This operation is equivalent to moving point E to point D. This is 
easily done if the curve Etot {N PW , Vi] is known (see Eq. © ). 

FIG. 4. Pressure in Si calculated at constant cut-off energy (E cu t=Q Ha with 2 special Appoints 
in the IBZ) for 51 lattice constants, (a) The open circles (o) represent the uncorrected values 
of pressure, whereas the solid diamonds (♦) illustrate the values corrected by the Froyen-Cohen 
technique, (b) The open diamonds (0) represent partially corrected values of pressure (see P* 
defined by Eq. (|3l"l) in the text), whereas the solid circles (•) are the values corrected by the S.H. 
technique. The graphs (a) and (b) point out that the S.H. correction technique (•) has two effects. 
The first is a shift of the uncorrected curve (o), also included in Froyen-Cohen technique (♦). 
Whereas the second is the cancellation of micro-curve jumps present in the partial correction ({>)• 
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FIG. 5. Total energy of Si calculated at constant cut-off energy (E cu t=6 Ha with 2 special 
fc-points in the IBZ) for 51 lattice constants. The open circles (o) represent the uncorrected values 
of total energy, whereas the solid circles (•) illustrate the values corrected by the S.H. technique. 
The solid curve and (resp.) the broken curve are obtained by a least-squares third-order fit to the 
corrected and (resp.) uncorrected data. The solid curve yields a lattice parameter of 10.23 Bohr, 
whereas the broken curve yields 10.17 Bohr. 

FIG. 6. Total energy of cubic BaTiO-3 calculated at constant cut-off energy (E cu t=20 Ha with 4 
special fc-points in the IBZ) for 4 lattice constants. The open circles (o) represent the uncorrected 
values of total energy, whereas the solid circles (•) are the values corrected by the S.H. technique. 
The corrected values of energy are joined by a least-squares third-order fit (solid curve). Without 
correction, the residual noise leads to incorrect values of the static equilibrium properties, as there 
are not enough data. The lattice parameter value is found to be 7.53 Bohr from the solid line. 

FIG. 7. Total energy and pressure of cubic BaTiC>3 calculated at constant cut-off energy 
(E cut =4:5 Ha with 10 special /c-points in the IBZ) for 8 lattice constants. The open circles (o) 
represent the uncorrected values (joined by a scattered dashed line), whereas the solid circles (•) 
illustrate the values corrected by the S.H. technique. The corrected values of energy and pressure 
are joined by a least-squares third-order fit (solid curves). For total energy, the values obtained 
with or without correction are nearly the same (a high energy cut-off has been used, which leads to 
well converged values of total energy). For pressure, instead, the values obtained with or without 
correction differ despite the high energy cut-off. The lattice parameter value is found to be 7.45 
Bohr from the solid lines. 
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FIG. 8. Semi-logarithmic plot of standard deviation \ °f the data as a function of the degree 
of the polynomial (obtained by least-squares fit). The standard deviation x, which is defined by 
Eq. (|3§| ) in the text, is a measure of the matching of a polynomial to the total energies of Si 
calculated at constant cut-off energy (E cu t=Q Ha with 2 special /c-points in the IBZ) for 51 lattice 
constants see Fig. §. The open circles (o) represent x calculated from uncorrected values of total 
energy, the open triangle (a) represent that calculated from the values corrected by Froyen-Cohen 
technique, whereas the solid circles (•) represent that calculated from the values corrected by the 
S.H. technique. The graph clearly illustrates the reduction of x with the correction. 

FIG. 9. Plot of standard deviations xe (top) and xp (bottom) of the energy and pressure 
data calculated at different scales at constant cut-off energy (E cut =6 Ha with 2 special /c-points 
in the IBZ) for 51 lattice constants. The white bars represent the standard deviation x calculated 
from uncorrected data, the gray bars illustrate the one calculated from the values corrected by 
Froyen-Cohen technique for pressure and by Francis-Payne technique for energy, whereas the black 
bars represent that calculated from the values corrected by the S.H. technique. The graph illus- 
trates that the reduction of % is more important for small scales than large ones. Note also that 
Froyen-Cohen technique does not reduce the noise in the pressure data. 
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FIG. 10. Total energy (top) and an (middle) of Si calculated at constant cut-off energy (E cu t=6 
Ha with 2 special fc-points in the IBZ) for 17 values of the ratio y = b/c in the case of an anistropic 
deformation which consists in compressing (or expanding) the unit cell along the [001] direction in 
order that the length of the three cubic directions are changed to a = b and c keeping the volume 
V unchanged. The open circles (o) represent the uncorrected values, whereas the solid circles (•) 
are the values corrected by the S.H. technique. The solid curves are obtained by least-squares 
third-order fit. The graph points out the effect of the correction, namely, the cancellation of the 
jumps between the micro-curves for energy, and a shift of the curve for stress. A decomposition 
of an into its different components is also presented (bottom). The open circles (o) represent the 
uncorrected values of isotropic part of an, whereas the solid circles (•) are the values corrected 
by the S.H. technique. The open diamonds (0) illustrate the values of the anisotropic part of an, 
that are not affected by the correction. 
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TABLES 

TABLE I. List of the different sets of lattice constants (expressed in atomic units) and the 
corresponding degree of the polynomial used to fit the data. 



Scale Points of the set Number of points Degree of the polynomial 

0.01 10.00 to 10.50 51 3 

0.025 10.00 to 10.50 21 3 

0.05 9.90 to 10.60 15 3 

0.1 9.50 to 11.00 16 3 

0.2 9.00 to 11.60 14 4 

0.5 7.00 to 13.00 13 8 
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TABLE II. Equilibrium lattice constant calculated at different scales (E cu t=6 Ha with 2 special 
/c-points in the IBZ) from the energy curve before (Before) and after correction with Francis-Payne 
technique (F-P) or scaling hypothesis technique (S.H.), and from the pressure curve before (Before) 
and after correction with Froyen-Cohen technique (F-C) or scaling hypothesis technique (S.H.). 
The values are expressed in Bohr, a is the standard deviation between the results obtained for the 
different scales. 





From energy 


From pressure 


Scale 


Before 


F-P 


S.H. 


Before 


F-C 


S.H. 


0.01 


10.1672 


10.2340 


10.2341 


10.0301 


10.2234 


10.2245 


0.025 


10.1678 


10.2349 


10.2351 


10.0312 


10.2236 


10.2243 


0.05 


10.1997 


10.2241 


10.2239 


10.0308 


10.2241 


10.2245 


0.1 


10.2136 


10.2266 


10.2262 


10.0295 


10.2236 


10.2224 


0.2 


10.2196 


10.2243 


10.2241 


10.0272 


10.2241 


10.2249 


0.5 


10.2287 


10.2188 


10.2188 


10.0280 


10.2256 


10.2241 


A* 


10.1994 


10.2271 


10.2270 


10.0294 


10.2241 


10.2241 


a 


2.6 10~ 2 


6.2 10~ 3 


6.3 10~ 3 


1.6 10~ 3 


8.1 10~ 4 


8.8 10~ 4 
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TABLE III. Equilibrium bulk modulus calculated at different scales (E cu t=6 Ha with 2 special 
fc-points in the IBZ) from the energy curve before (Before) and after correction with Francis-Payne 
technique (F-P) or scaling hypothesis technique (S.H.), and from the pressure curve before (Before) 
and after correction with Froyen-Cohen technique (F-C) or scaling hypothesis technique (S.H.). 
The values are expressed in Mbar. a is the standard deviation between the results obtained for 
the different scales. 





From energy 


From pressure 


Scale 


Before 


F-P 


S.H. 


Before 


F-C 


S.H. 


0.01 


1.181 


0.8243 


0.8233 


1.071 


0.9849 


0.9419 


0.025 


1.079 


0.8390 


0.8368 


1.088 


0.9823 


0.9417 


0.05 


0.8326 


0.8860 


0.8859 


1.123 


0.9560 


0.9426 


0.1 


1.0030 


0.9709 


0.9710 


1.127 


0.9482 


0.9412 


0.2 


0.9546 


0.9630 


0.9601 


1.115 


0.9333 


0.9285 


0.5 


0.9539 


0.9485 


0.9481 


1.108 


0.9321 


0.9411 


A* 


1.0007 


0.9053 


0.9042 


1.1053 


0.9561 


0.9395 


a 


1.2 1Q- 1 


6.5 10~ 2 


6.5 10~ 2 


2.2 10~ 2 


2.3 10~ 2 


5.4 10~ 3 
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TABLE IV. Equilibrium lattice constant calculated for the different micro-curves (E cu t=6 Ha 
with 2 special fe-points in the IBZ) from the energy curve before (Before) and after correction 
with Francis-Payne technique (F-P) or scaling hypothesis technique (S.H.), and from the pressure 
curve before (Before) and after correction with Froyen-Cohen technique (F-C) or scaling hypothesis 
technique (S.H.). The values are expressed in Bohr, a is the standard deviation between the results 
obtained for the different scales. 





From energy 


From pressure 


TVTi pro- n i rvp 


JJC1UI c 


F-P 

_L — 1 


D .11. 


JJC1U1 c 


F-C 

1 V_y 


S H 

kJ .11. 


1 


10.0263 


10.2000 


10.2287 


10.0263 


10.1996 


10.2279 


2 


10.0380 


10.2138 


10.2262 


10.0380 


10.2135 


10.2258 


3 


10.0423 


10.2188 


10.2260 


10.0425 


10.2188 


10.2260 


4 


10.0537 


10.2247 


10.2241 


10.0473 


10.2245 


10.2241 


5 


10.0529 


10.2328 


10.2224 


10.0544 


10.2328 


10.2226 


6 


10.0577 


10.2400 


10.2216 


10.0605 


10.2404 


10.2222 


7 


10.0671 


10.2527 


10.2200 


10.0711 


10.2534 


10.2213 


A* 


10.0482 


10.2261 


10.2241 


10.0485 


10.2261 


10.2241 


a 


1.4 10~ 2 


1.7 10" 2 


3.1 10~ 3 


1.5 10~ 2 


1.8 10~ 2 


2.4 10~ 3 
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TABLE V. Equilibrium bulk modulus calculated for the different micro-curves (E cu t=6 Ha 
with 2 special fe-points in the IBZ) from the energy curve before (Before) and after correction 
with Francis-Payne technique (F-P) or scaling hypothesis technique (S.H.), and from the pressure 
curve before (Before) and after correction with Froyen-Cohen technique (F-C) or scaling hypothesis 
technique (S.H.). The values are expressed in Mbar. a is the standard deviation between the results 
obtained for the differentscales. 





From energy 


From pressure 


Mi cro-fi i rve 




F-P 


S H 




F-C 

X v_ 


S H 

kJ .XX . 


1 


1.244 


1.080 


0.9178 


1.244 


1.087 


0.9294 


2 


1.226 


1.063 


0.9203 


1.227 


1.069 


0.9278 


3 


1.215 


1.062 


0.9253 


1.220 


1.062 


0.9261 


4 


1.278 


1.071 


0.9428 


1.211 


1.056 


0.9281 


5 


1.125 


1.044 


0.9260 


1.201 


1.047 


0.9295 


6 


1.165 


1.031 


0.9201 


1.189 


1.037 


0.9285 


7 


1.142 


1.012 


0.9143 


1.171 


1.022 


0.9276 




1.1993 


1.0518 


0.9238 


1.2090 


1.0543 


0.9281 


a 


5 1(T 2 


2 10~ 2 


i io- 2 


2 10~ 2 


2 10~ 2 


1 10~ 3 
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TABLE VI. Equilibrium lattice constant calculated with different cut-off energy at scale 0.1 
(with 2 special /c-points in the IBZ) from the energy curve before (Before) and after correction with 
Francis-Payne technique (F-P) or scaling hypothesis technique (S.H.), and from the pressure curve 
before (Before) and after correction with Froyen-Cohen technique (F-C) or scaling hypothesis 
technique (S.H.). The values are expressed in Bohr. The values between brackets indicate the 
relative error with respect to the fully converged value at 15 Ha, for the same technique. 





From energy 


From pressure 


E C ut 


Before 


F-P 


S.H. 


Before 


F-C 


S.H. 


3 


10.1121 


10.0935 


10.0935 


9.6899 


10.0856 


10.0854 




(0.8%) 


(0.9%) 


(0.9%) 


(5%) 


(1%) 


(1%) 


6 


10.2135 


10.2266 


10.2262 


10.0295 


10.2236 


10.224 




(0.2%) 


(0.4%) 


(0.4%) 


(2%) 


(0.4%) 


(0.4%) 


10 


10.1897 


10.1894 


10.1894 


10.1771 


10.1867 


10.1869 




(<0.1%) 


(«0.1%) 


(<0.1%) 


(«0.1%) 


(<0.1%) 


(«0.1%) 


15 


10.1886 


10.1888 


10.1888 


10.1716 


10.1856 


10.1856 
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TABLE VII. Equilibrium bulk modulus calculated with different cut-off energy at scale 0.1 
(with 2 special /c-points in the IBZ) from the energy curve before (Before) and after correction 
with Francis-Payne technique (F-P) or scaling hypothesis technique (S.H.), and from the pressure 
curve before (Before) and after correction with Froyen-Cohen technique (F-C) or scaling hypothesis 
technique (S.H.). The values are expressed in Mbar. The values between brackets indicate the 
relative error with respect to the fully converged value at 15 Ha, for the same technique. 





From energy 


From pressure 


E C ut 


Before 


F-P 


S.H. 


Before 


F-C 


S.H. 


3 


1.126 


1.169 


1.169 


1.620 


1.154 


1.154 




(14%) 


(18%) 


(18%) 


(64%) 


(19%) 


(19%) 


6 


1.003 


0.9709 


0.9710 


1.1270 


0.9482 


0.9412 




(2%) 


(2%) 


(2%) 


(15%) 


(2%) 


(3%) 


10 


0.9889 


0.9902 


0.9902 


0.9805 


0.9712 


0.9714 




(0.2%) 


(0.3%) 


(0.3%) 


(0.2%) 


(0.2%) 


(0.2%) 


15 


0.9875 


0.9872 


0.9872 


0.9826 


0.9692 


0.9692 
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10.0 10.1 10.2 10.3 10.4 10.5 

a (Bohr) 




10.0 10.1 10.2 10.3 10.4 10.5 

a (Bohr) 




10.0 10.1 10.2 10.3 10.4 10.5 

a (Bohr) 



